Model & visualize parameter estimates in R

ROIs from the Craddock et al. (2012) parcellation atlas

Mean parameter estimates were extracted from parcel 292 and 116. The parcellation atlas can be found in `functional-workshop/data/ROIs/craddock_all.nii.gz. This atlas has multiple volumes within the 4D file, and we extracted from the K=400 atlas, which is the 31st volume in AFNI (index = 0 in AFNI).

Extract mean parameter estimates

Run bash script to calculate mean parameter estimates for each subject, wave, and condition contrast (condition > rest) within each ROI using AFNI 3dmaskave.

Path to bash script: functional-workshop/code/ROI_analysis/extract_parameterEstimates.sh

#!/bin/bash
. ~/.bashrc

# This script extracts mean parameter estimates and SDs within an ROI or parcel
# from subject FX condition contrasts (condition > rest) for each wave. Output is 
# saved as a text file in the output directory.

# Set paths and variables
# ------------------------------------------------------------------------------------------
# paths

con_dir='/Volumes/psych-cog/dsnlab/MDC/functional-workshop/data/FX_models' #fx contrast directory
atlas_dir='/Volumes/psych-cog/dsnlab/MDC/functional-workshop/data/ROIs' #roi/parcellation atlas directory 
output_dir='/Volumes/psych-cog/dsnlab/MDC/functional-workshop/results/ROI_analysis' #roi/parcellation output directory
rx_model='/Volumes/psych-cog/dsnlab/MDC/functional-workshop/results/AFNI/all+tlrc' #rx model (for atlas alignment only)

# variables
subjects=`cat /Volumes/psych-cog/dsnlab/MDC/functional-workshop/data/subject_list.txt`
parcellation_atlas=(craddock_all.nii.gz) #roi/parcellation atlas file
parcellation_map=(31) #parcellation map number (if applicable)
aligned_parcellation_map=(aligned_craddock_400) #aligned roi/parcellation map name
aligned_parcellation_num=(116 292) #parcellation number(s) to extract from; use $(seq 1 N) where N is the total number of parcels to extract from all
waves=(t1 t2 t3) #waves or task names
fx_cons=(con_0001 con_0002 con_0003 con_0004) #fx con files to extract from

# Align roi/parcellation map to data
# ------------------------------------------------------------------------------------------
echo "aligning parcellation map"
if [ -f $atlas_dir/${aligned_parcellation_map}+tlrc.BRIK ]; then
    echo "aligned parcellation map already exists"
else 
3dAllineate -source $atlas_dir/$parcellation_atlas[$parcellation_map] -master $rx_model -final NN -1Dparam_apply '1D: 12@0'\' -prefix $atlas_dir/$aligned_parcellation_map
fi

# Extract mean parameter estimates and SDs for each subject, wave, contrast, and roi/parcel
# ------------------------------------------------------------------------------------------

for sub in ${subjects[@]}; do 
    for wave in ${waves[@]}; do 
        for con in ${fx_cons[@]}; do 
            for num in ${aligned_parcellation_num[@]}; do 
                echo ${sub} ${wave} ${con} ${num} `3dmaskave -sigma -quiet -mrange $num $num -mask $atlas_dir/${aligned_parcellation_map}+tlrc $con_dir/${sub}_${wave}_${con}.nii` >> $output_dir/parameterEstimates.txt
            done
        done
    done
done

The output will be saved in a text file functional-workshop/results/ROI_analysis/parameterEstimates.txt

library(dplyr)
library(knitr)

read.table('/Volumes/psych-cog/dsnlab/MDC/functional-workshop/results/ROI_analysis/parameterEstimates.txt', sep = "", fill = TRUE, stringsAsFactors=FALSE) %>%
  head(10) %>%
  kable(format = 'pandoc')

V1	V2	V3	V4	V5	V6
s001	t1	con_0001	116	-0.5871270	0.274073
s001	t1	con_0001	292	-0.3261990	0.331749
s001	t1	con_0002	116	-0.0216224	0.342096
s001	t1	con_0002	292	-0.0078532	0.387301
s001	t1	con_0003	116	-0.5910870	0.330697
s001	t1	con_0003	292	-0.2818540	0.396804
s001	t1	con_0004	116	-0.0967547	0.305696
s001	t1	con_0004	292	-0.0509896	0.374470
s001	t2	con_0001	116	-0.6617310	0.311137
s001	t2	con_0001	292	-0.4697610	0.275190

Load packages

library(tidyverse)
library(lme4)
library(lmerTest)
library(wesanderson)

Load data

# load parameter estimate .csv file
data = read.table('/Volumes/psych-cog/dsnlab/MDC/functional-workshop/results/ROI_analysis/parameterEstimates.txt', sep = "", fill = TRUE, stringsAsFactors=FALSE)

# load age covariates and rename variables
age = read.csv('/Volumes/psych-cog/dsnlab/MDC/functional-workshop/data/covariates/age.csv') %>%
  rename("subjectID" = Subj,
         "wave" = wavenum)

Tidy data

Specify your variables names and levels

# tidy raw data
data1 = data %>% 
  # rename variables
  rename('subjectID' = V1,
         'wave' = V2,
         'con' = V3,
         'parcellation' = V4,
         'beta' = V5,
         'sd' = V6) %>%
  # convert con file names to condition names
  mutate(target = ifelse(con %in% c('con_0001', 'con_0002'), 'self', 'other'), 
         domain = ifelse(con %in% c('con_0001', 'con_0003'), 'academic', 'social'), 
  # change data type to factor
         parcellation = as.factor(parcellation),
         target = as.factor(target),
         domain = as.factor(domain)) %>%
  # change to integer
  extract(wave, 'wave', 't([0-3]{1})') %>%
  mutate(wave = as.integer(wave))

Merge data, add age to the data frame and center

merged = left_join(data1, age, by = c('subjectID', 'wave')) %>%
  mutate(age_c = age-mean(age, na.rm=TRUE))

# print data frame header
merged %>%
  head(16) %>%
  kable(format = 'pandoc')

subjectID	wave	con	parcellation	beta	sd	target	domain	age	age_c
s001	1	con_0001	116	-0.5871270	0.274073	self	academic	10.41530	-2.0912751
s001	1	con_0001	292	-0.3261990	0.331749	self	academic	10.41530	-2.0912751
s001	1	con_0002	116	-0.0216224	0.342096	self	social	10.41530	-2.0912751
s001	1	con_0002	292	-0.0078532	0.387301	self	social	10.41530	-2.0912751
s001	1	con_0003	116	-0.5910870	0.330697	other	academic	10.41530	-2.0912751
s001	1	con_0003	292	-0.2818540	0.396804	other	academic	10.41530	-2.0912751
s001	1	con_0004	116	-0.0967547	0.305696	other	social	10.41530	-2.0912751
s001	1	con_0004	292	-0.0509896	0.374470	other	social	10.41530	-2.0912751
s001	2	con_0001	116	-0.6617310	0.311137	self	academic	13.48767	0.9810956
s001	2	con_0001	292	-0.4697610	0.275190	self	academic	13.48767	0.9810956
s001	2	con_0002	116	-0.3545030	0.338036	self	social	13.48767	0.9810956
s001	2	con_0002	292	-0.1359760	0.424825	self	social	13.48767	0.9810956
s001	2	con_0003	116	-0.5791710	0.404919	other	academic	13.48767	0.9810956
s001	2	con_0003	292	-0.5360000	0.272914	other	academic	13.48767	0.9810956
s001	2	con_0004	116	-0.2961900	0.273353	other	social	13.48767	0.9810956
s001	2	con_0004	292	-0.2756520	0.265749	other	social	13.48767	0.9810956

Remove missing data to run LME models

data.complete = merged %>%
  na.omit(.)

# print number of rows
cat('raw data: ', nrow(merged))

## raw data:  1944

cat('complete data: ', nrow(data.complete))

## complete data:  1296

Run LME models within parcel 292 and compare

Predict parameter estimates from task conditions (target and domain) and age within parcel 292.

Linear effect of age, random intercepts only

model.1 = lmer(beta ~ target*domain*age_c + (1 | subjectID), data=filter(data.complete, parcellation == 292))
summary(model.1)

## Linear mixed model fit by REML t-tests use Satterthwaite approximations
##   to degrees of freedom [lmerMod]
## Formula: beta ~ target * domain * age_c + (1 | subjectID)
##    Data: filter(data.complete, parcellation == 292)
## 
## REML criterion at convergence: 1106.1
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.0561 -0.4861  0.0324  0.4626  4.3648 
## 
## Random effects:
##  Groups    Name        Variance Std.Dev.
##  subjectID (Intercept) 0.1008   0.3175  
##  Residual              0.2582   0.5082  
## Number of obs: 648, groups:  subjectID, 81
## 
## Fixed effects:
##                                Estimate Std. Error        df t value
## (Intercept)                    -0.07709    0.05422 199.40000  -1.422
## targetself                      0.03147    0.05663 554.60000   0.556
## domainsocial                    0.09379    0.05663 554.60000   1.656
## age_c                          -0.03904    0.01658 584.80000  -2.354
## targetself:domainsocial         0.09714    0.08008 554.60000   1.213
## targetself:age_c                0.06418    0.02278 554.60000   2.817
## domainsocial:age_c              0.02990    0.02278 554.60000   1.312
## targetself:domainsocial:age_c  -0.06094    0.03222 554.60000  -1.891
##                               Pr(>|t|)   
## (Intercept)                    0.15660   
## targetself                     0.57861   
## domainsocial                   0.09824 . 
## age_c                          0.01889 * 
## targetself:domainsocial        0.22563   
## targetself:age_c               0.00502 **
## domainsocial:age_c             0.19000   
## targetself:domainsocial:age_c  0.05911 . 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) trgtsl dmnscl age_c  trgts: trgt:_ dmns:_
## targetself  -0.522                                          
## domainsocil -0.522  0.500                                   
## age_c        0.092 -0.052 -0.052                            
## trgtslf:dmn  0.369 -0.707 -0.707  0.037                     
## trgtslf:g_c -0.040  0.076  0.038 -0.687 -0.054              
## domnscl:g_c -0.040  0.038  0.076 -0.687 -0.054  0.500       
## trgtslf:d:_  0.028 -0.054 -0.054  0.486  0.076 -0.707 -0.707

Linear effect of age, random intercepts and age slopes

model.2 = lmer(beta ~ target*domain*age_c + (1 + age_c | subjectID), data=filter(data.complete, parcellation == 292))
summary(model.2)

## Linear mixed model fit by REML t-tests use Satterthwaite approximations
##   to degrees of freedom [lmerMod]
## Formula: beta ~ target * domain * age_c + (1 + age_c | subjectID)
##    Data: filter(data.complete, parcellation == 292)
## 
## REML criterion at convergence: 1101.3
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.8144 -0.4814  0.0320  0.4554  4.3703 
## 
## Random effects:
##  Groups    Name        Variance Std.Dev. Corr 
##  subjectID (Intercept) 0.092839 0.30469       
##            age_c       0.001157 0.03402  -0.60
##  Residual              0.252104 0.50210       
## Number of obs: 648, groups:  subjectID, 81
## 
## Fixed effects:
##                                Estimate Std. Error        df t value
## (Intercept)                    -0.07347    0.05281 204.40000  -1.391
## targetself                      0.03147    0.05595 516.30000   0.562
## domainsocial                    0.09379    0.05595 516.30000   1.676
## age_c                          -0.03807    0.01680 307.30000  -2.266
## targetself:domainsocial         0.09714    0.07913 516.30000   1.228
## targetself:age_c                0.06418    0.02251 516.30000   2.851
## domainsocial:age_c              0.02990    0.02251 516.30000   1.328
## targetself:domainsocial:age_c  -0.06094    0.03184 516.30000  -1.914
##                               Pr(>|t|)   
## (Intercept)                    0.16569   
## targetself                     0.57405   
## domainsocial                   0.09429 . 
## age_c                          0.02412 * 
## targetself:domainsocial        0.22013   
## targetself:age_c               0.00453 **
## domainsocial:age_c             0.18475   
## targetself:domainsocial:age_c  0.05617 . 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) trgtsl dmnscl age_c  trgts: trgt:_ dmns:_
## targetself  -0.530                                          
## domainsocil -0.530  0.500                                   
## age_c       -0.003 -0.051 -0.051                            
## trgtslf:dmn  0.375 -0.707 -0.707  0.036                     
## trgtslf:g_c -0.040  0.076  0.038 -0.670 -0.054              
## domnscl:g_c -0.040  0.038  0.076 -0.670 -0.054  0.500       
## trgtslf:d:_  0.028 -0.054 -0.054  0.474  0.076 -0.707 -0.707

Compare models

model.1: beta ~ target * domain * age_c + (1 | subjectID)

model.2: beta ~ target * domain * age_c + (1 + age_c | subjectID)

anova(model.1, model.2) %>%
  `row.names<-`(c('model.1', 'model.2')) %>%
  kable()

	Df	AIC	BIC	logLik	deviance	Chisq	Chi Df	Pr(>Chisq)
model.1	10	1083.649	1128.388	-531.8247	1063.649	NA	NA	NA
model.2	12	1082.726	1136.413	-529.3631	1058.726	4.923206	2	0.0852981

Adding age as a random effect does not significantly improve the model fit.

Visualize raw data

# set color palette
palette = wes_palette("Zissou", 2, type = "continuous")

Plot fitted curves for parcels 292 and 116

Main effect of target

ggplot(data.complete, aes(x = age, 
                          y = beta, 
                          group = interaction(subjectID, target, domain), 
                          color = target)) +
  geom_point(size = .5, alpha = .1) + 
  geom_line(alpha = .1) + 
  geom_line(aes(group=target), size = 1.5, stat = 'smooth', method = 'lm', formula = y ~ poly(x,2)) + 
  facet_wrap(~parcellation, ncol = 2) +
  geom_hline(yintercept = 0, color = 'gray') +
  scale_color_manual(breaks = c('self', 'other'), values = c(self=palette[2], other=palette[1])) +
  scale_x_continuous(breaks=c(10,13,16)) +
  coord_cartesian(ylim=c(-1,1)) +
  theme_minimal(base_size = 18)

Interaction between target and domain

ggplot(data.complete, aes(x = age, 
                          y = beta, 
                          group = interaction(subjectID, target, domain), 
                          color = target, 
                          linetype = domain)) +
  geom_point(size = .5, alpha = .1) + 
  geom_line(alpha = .1) + 
  geom_line(aes(group=interaction(target,domain)), size = 1.5, stat = 'smooth', method = 'lm', formula = y ~ poly(x,2)) + 
  facet_wrap(~parcellation, ncol = 2) +
  geom_hline(yintercept = 0, color = 'gray')+
  scale_color_manual(breaks = c('self', 'other'), values = c(self=palette[2], other=palette[1]))+
  scale_x_continuous(breaks=c(10,13,16)) +
  coord_cartesian(ylim=c(-1,1)) +
  theme_minimal(base_size = 18)

Plot LOESS curves for parcels 292 and 116

Main effect of target

ggplot(data.complete, aes(x = age, 
                          y = beta, 
                          group = interaction(subjectID, target, domain), 
                          color = target, 
                          linetype = domain)) +
  geom_point(size = .5, alpha = .1) + 
  geom_line(alpha = .1) + 
  geom_line(aes(group=target), size = 1.5, stat = 'smooth', method = 'loess') + 
  facet_wrap(~parcellation, ncol = 2) +
  geom_hline(yintercept = 0, color = 'gray')+
  scale_color_manual(breaks = c('self', 'other'), values = c(self=palette[2], other=palette[1]))+
  scale_x_continuous(breaks=c(10,13,16)) +
  coord_cartesian(ylim=c(-1,1)) +
  theme_minimal(base_size = 18)

Interaction between target and domain

ggplot(data.complete, aes(x = age, 
                          y = beta, 
                          group = interaction(subjectID, target, domain), 
                          color = target, 
                          linetype = domain)) +
  geom_point(size = .5, alpha = .1) + 
  geom_line(alpha = .1) + 
  geom_line(aes(group=interaction(target,domain)), size = 1.5, stat = 'smooth', method = 'loess') + 
  facet_wrap(~parcellation, ncol = 2) +
  geom_hline(yintercept = 0, color = 'gray')+
  scale_color_manual(breaks = c('self', 'other'), values = c(self=palette[2], other=palette[1]))+
  scale_x_continuous(breaks=c(10,13,16)) +
  coord_cartesian(ylim=c(-1,1)) +
  theme_minimal(base_size = 18)

Visualize predicted values from model.1

Linear effect of age, random intercepts only

model.1: beta ~ target * domain * age_c + (1 | subjectID)

Plot fitted curves for parcels 292 and 116

# extract random effects formula from model.2 and reconstruct it to use with the `predict` function
REFormulaString = as.character(findbars(model.1@call$formula)[[1]])
REFormula = as.formula(paste0('~(', REFormulaString[[2]], REFormulaString[[1]], REFormulaString[[3]], ')'))

# get expected values for each observation based on model.2
data.complete$expected.1 <- predict(model.1, newdata = data.complete, re.form=REFormula)
data.complete$expected_mean.1 <- predict(model.1, newdata = data.complete, re.form=NA)

Main effect of target

ggplot(data.complete, aes(x = age, 
                          y = expected.1, 
                          group = interaction(subjectID, target, domain), 
                          color = target)) +
  geom_point(size = .5, alpha = .1) + 
  geom_line(alpha = .1) + 
  geom_line(aes(y = expected_mean.1, group=target), size = 1.5, stat = 'smooth', method = 'lm', formula = y ~ poly(x,2)) + 
  facet_wrap(~parcellation, ncol = 2) +
  geom_hline(yintercept = 0, color = 'gray') +
  scale_color_manual(breaks = c('self', 'other'), values = c(self=palette[2], other=palette[1])) +
  scale_x_continuous(breaks=c(10,13,16)) +
  coord_cartesian(ylim=c(-1,1)) +
  theme_minimal(base_size = 18)

Interaction between target and domain

ggplot(data.complete, aes(x = age, 
                          y = expected.1, 
                          group = interaction(subjectID, target, domain), 
                          color = target, 
                          linetype = domain)) +
  geom_point(size = .5, alpha = .1) + 
  geom_line(alpha = .1) + 
  geom_line(aes(y = expected_mean.1, group=interaction(target,domain)), size = 1.5, stat = 'smooth', method = 'lm', formula = y ~ poly(x,2)) + 
  facet_wrap(~parcellation, ncol = 2) +
  geom_hline(yintercept = 0, color = 'gray')+
  scale_color_manual(breaks = c('self', 'other'), values = c(self=palette[2], other=palette[1]))+
  scale_x_continuous(breaks=c(10,13,16)) +
  coord_cartesian(ylim=c(-1,1)) +
  theme_minimal(base_size = 18)

Visualize predicted values from model.2

Linear effect of age, random intercepts and age slopes

model.2: beta ~ target * domain * age_c + (1 + age_c | subjectID)

Plot fitted curves for parcels 292 and 116

# extract random effects formula from model.2 and reconstruct it to use with the `predict` function
REFormulaString = as.character(findbars(model.2@call$formula)[[1]])
REFormula = as.formula(paste0('~(', REFormulaString[[2]], REFormulaString[[1]], REFormulaString[[3]], ')'))

# get expected values for each observation based on model.2
data.complete$expected.2 <- predict(model.2, newdata = data.complete, re.form=REFormula)
data.complete$expected_mean.2 <- predict(model.2, newdata = data.complete, re.form=NA)

Main effect of target

ggplot(data.complete, aes(x = age, 
                          y = expected.2, 
                          group = interaction(subjectID, target, domain), 
                          color = target)) +
  geom_point(size = .5, alpha = .1) + 
  geom_line(alpha = .1) + 
  geom_line(aes(y = expected_mean.2, group=target), size = 1.5, stat = 'smooth', method = 'lm', formula = y ~ poly(x,2)) + 
  facet_wrap(~parcellation, ncol = 2) +
  geom_hline(yintercept = 0, color = 'gray') +
  scale_color_manual(breaks = c('self', 'other'), values = c(self=palette[2], other=palette[1])) +
  scale_x_continuous(breaks=c(10,13,16)) +
  coord_cartesian(ylim=c(-1,1)) +
  theme_minimal(base_size = 18)

Interaction between target and domain

ggplot(data.complete, aes(x = age, 
                          y = expected.2, 
                          group = interaction(subjectID, target, domain), 
                          color = target, 
                          linetype = domain)) +
  geom_point(size = .5, alpha = .1) + 
  geom_line(alpha = .1) + 
  geom_line(aes(y = expected_mean.2, group=interaction(target,domain)), size = 1.5, stat = 'smooth', method = 'lm', formula = y ~ poly(x,2)) + 
  facet_wrap(~parcellation, ncol = 2) +
  geom_hline(yintercept = 0, color = 'gray')+
  scale_color_manual(breaks = c('self', 'other'), values = c(self=palette[2], other=palette[1]))+
  scale_x_continuous(breaks=c(10,13,16)) +
  coord_cartesian(ylim=c(-1,1)) +
  theme_minimal(base_size = 18)